Comparisons of various approaches to determine optimal placement of
restoration areas given an underlying habitat shp file
using Lizard Island example from Allen Coral Atlas.
shp files, filter and remove polygons less than
target plot hectare size (1ha).Red inset box = subset map used below (note - spatial data in projected coordinate system 20353)
library(tidyverse)
library(sf)
library(future.apply)
library(tictoc)
library(lwgeom)
library(tmap)
set.seed(123)
# Import reef polygons, create union (or returns multiple intersections
geomorphic_seeding <- st_read("/Users/rof011/spatialtools/apps/remove-zones/www/geomorphic.geojson") %>%
mutate(area=as.numeric(st_area(.))) |>
st_transform(20353) |>
filter(area > 10000) |>
mutate(n=round(as.numeric(area)/1000)) |>
st_transform(4326) |>
st_crop(st_bbox(c(xmin = 145.42, xmax = 145.48, ymin = -14.72, ymax = -14.64))) |>
st_transform(20353)
## Reading layer `Lizard_Eyrie' from data source
## `/Users/rof011/spatialtools/apps/remove-zones/www/geomorphic.geojson'
## using driver `GeoJSON'
## Simple feature collection with 736 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 145.3519 ymin: -14.73534 xmax: 145.4974 ymax: -14.64425
## Geodetic CRS: WGS 84
geomorphic_seeding_filtered <- geomorphic_seeding %>%
filter(class %in% c("Outer Reef Flat", "Reef Slope", "Back Reef Slope"))
inset <- st_bbox(c(xmin = 145.44, xmax = 145.455, ymin = -14.697, ymax = -14.692)) |>
st_bbox()
inset_sf <- inset |> st_as_sfc() |> st_set_crs(4326)
habitat_pal <- c("Plateau" = "cornsilk2", "Back Reef Slope" = "darkcyan",
"Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
"Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
"Reef Crest" = "coral3", "Shallow Lagoon" = "turquoise3",
"Deep Lagoon" = "turquoise4")
habitat_pal_filtered <- c("Back Reef Slope" = "darkcyan",
"Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
"Outer Reef Flat" = "darkgoldenrod2")
tmap_mode("plot")
tm_basemap("Esri.WorldImagery", alpha=0.6) +
tm_shape(geomorphic_seeding) +
tm_polygons(fill="class",
lwd=0,
fill_alpha=0.7,
fill.scale=tm_scale_categorical(values=habitat_pal)) +
tm_shape(inset_sf) +
tm_polygons(fill=NA,
col="red",
lwd=2)
Determine optimum plot 1ha locations using a regular parallel grid (100m x 100m).
geomorphic_seeding_union <- geomorphic_seeding %>%
filter(class %in% c("Outer Reef Flat", "Reef Slope", "Back Reef Slope", "Sheltered Reef Slope")) %>%
st_snap_to_grid(size = 0.01) %>%
st_make_valid() %>%
st_buffer(dist=0.1) %>%
group_by() %>%
summarise(area = sum(area)) %>%
ungroup()
geomorphic_seeding_buffered <- st_buffer(geomorphic_seeding_union, dist = 100)
bbox <- st_bbox(geomorphic_seeding_buffered)
grid <- st_make_grid(bbox, cellsize = c(100, 100), square = TRUE)
grid_sf <- st_sf(geometry = grid)
contains_properly <- st_contains_properly(geomorphic_seeding_union, grid_sf)
indices <- unlist(contains_properly)
indices_numeric <- as.numeric(indices)
grid_sf_whole <- grid_sf[indices_numeric, ]
# Visualization with tmap
tm_basemap("Esri.WorldImagery") +
tm_shape(geomorphic_seeding_filtered, bbox=inset) +
tm_polygons(fill = "class",
fill.legend = tm_legend(title="Habitat", position = tm_pos_out("right", "center")),
fill.scale=tm_scale_categorical(values=habitat_pal_filtered),
lwd = 0.2, col="black", fill_alpha = 1) +
tm_shape(grid_sf, bbox = inset, crs = 20353) +
tm_polygons(fill = NA, col = "darkred", lwd = 0.5) +
tm_shape(grid_sf_whole, bbox = inset, crs = 20353) +
tm_polygons(fill = NA, col = "white", lwd = 1) +
tm_title("1) Grid polygons", size=1)
to iterate, use lapply to shift the grid by 10m
increments in lon and lat (or finer
resolution):
library(sf)
library(dplyr)
library(tmap)
# Ensure CRS is set to 20353
crs_value <- 20353
# Define the steps for shifting the grid
step_size <- 10
shifts <- seq(0, 90, by = step_size) # 10m steps up to 90m
bbox <- st_bbox(geomorphic_seeding_union)
# Function to shift grid and filter based on st_contains_properly
generate_shifted_grid <- function(shift_x, shift_y) {
# Shift the bounding box
shifted_bbox <- bbox
shifted_bbox['xmin'] <- bbox['xmin'] + shift_x
shifted_bbox['xmax'] <- bbox['xmax'] + shift_x
shifted_bbox['ymin'] <- bbox['ymin'] + shift_y
shifted_bbox['ymax'] <- bbox['ymax'] + shift_y
# Ensure that the shifted_bbox values are valid
if (any(is.na(shifted_bbox))) {
stop("Shifted bounding box contains NA values, check your shifts.")
}
# Generate the grid and set the CRS
grid <- st_make_grid(shifted_bbox, cellsize = c(100, 100), square = TRUE)
grid_sf <- st_sf(geometry = grid, crs = crs_value)
# Check if grid cells are fully within the geomorphic_seeding_union
contains_properly <- st_contains_properly(geomorphic_seeding_union, grid_sf)
indices <- unlist(contains_properly)
# Add a color column based on whether the grid cells are fully contained
grid_sf <- grid_sf %>%
mutate(color = ifelse(row_number() %in% indices, "white", "darkred"))
return(grid_sf)
}
# Apply the function over all combinations of shifts in x and y directions
grids_list <- lapply(shifts, function(x_shift) {
lapply(shifts, function(y_shift) {
generate_shifted_grid(x_shift, y_shift)
})
})
# Flatten the nested list of sf objects into a single list
flattened_list <- do.call(c, grids_list)
# Combine all sf objects into a single sf object
grids_list <- do.call(rbind, flattened_list)
gridplots <- tm_basemap("Esri.WorldImagery") +
tm_shape(geomorphic_seeding_filtered, bbox=inset) +
tm_polygons(fill = "class",
fill.legend = tm_legend(title="Habitat", position = tm_pos_out("right", "center")),
fill.scale=tm_scale_categorical(values=habitat_pal_filtered),
lwd = 0.2, col="black", fill_alpha = 1) +
tm_shape(grids_list, bbox = inset, crs = 20353) +
tm_polygons(fill = NA, col = "color", lwd = "color",
col.legend = tm_legend(title="Habitat", position = tm_pos_out("right", "center")),
lwd.scale=tm_scale_categorical(values=c(0, 1))) +
tm_title("1) Gridded plots")
gridplots
Determine optimum plot 1ha locations using polygons with randomly seeded centroids, keep polygons that are 100% overlap with habitat type.
The number of seed points (centroids for polygons) is determined by the area of each polygon (i.e. larger polygons = more points)
set.seed(101)
#remotes::install_github("r-tmap/tmap.deckgl")
# lapply sampling process
seeded_points <- st_sf(do.call(rbind, lapply(seq_len(nrow(geomorphic_seeding_filtered)), function(i) {
polygon <- geomorphic_seeding_filtered[i, ]
points <- st_sample(polygon, size = polygon$n, type = "random")
st_sf(class = polygon$class, geometry = points)
}))) |> st_transform(20353)
tm_basemap("Esri.WorldImagery") +
tm_shape(geomorphic_seeding_filtered) +
tm_polygons(fill = "class",
fill.legend = tm_legend(title="Habitat", position = tm_pos_out("right", "center")),
fill.scale=tm_scale_categorical(values=habitat_pal_filtered),
lwd = 0.2, col="black", fill_alpha = 1) +
tm_shape(seeded_points) +
tm_dots(size=0.01,
fill="white",
fill_alpha=0.5) +
tm_shape(inset_sf) +
tm_polygons(fill=NA,
col="red",
lwd=2) +
tm_title("Random seed centroids for polygons")
Generate 100 x 100m polygons around each point:
# Define the width and length of the rectangle
width <- 100 # Example width (in the same units as your CRS)
length <- 100 # Example length (in the same units as your CRS)
# Use lapply to create the buffered polygons
buffered_polygons <- lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
return(polygon)
})
# Combine the results into a single sf object
buffered_polygons_sf <- sf::st_sfc(buffered_polygons, crs = sf::st_crs(seeded_points)) |> st_as_sf()
tm_basemap("Esri.WorldImagery") +
tm_shape(geomorphic_seeding_filtered, bbox=inset) +
tm_polygons(fill = "class",
fill.legend = tm_legend(title="Habitat", position = tm_pos_out("right", "center")),
fill.scale=tm_scale_categorical(values=habitat_pal_filtered),
lwd = 0.2, col="black", fill_alpha = 1) +
tm_shape(buffered_polygons_sf,
bbox=inset,
crs=20353) +
tm_polygons(fill=NA,
col="white",
lwd=0.8) +
tm_title("Random parallel polygons", size=1)
Make a function to apply plots across points using
lapply and retain only polygons with 100% overlap with reef
habitats as with gridded approach:
library(tictoc)
tic()
# lapply sampling process
seeded_points <- st_sf(do.call(rbind, lapply(seq_len(nrow(geomorphic_seeding_filtered)), function(i) {
polygon <- geomorphic_seeding_filtered[i, ]
points <- st_sample(polygon, size = polygon$n, type = "random")
st_sf(class = polygon$class, geometry = points)
}))) |> st_transform(20353)
# Define the width and length of the rectangle
width <- 100
length <- 100
# Use lapply to create the buffered polygons
buffered_polygons_filtered <- lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
return(polygon)
})
# Combine the results into a single sf object
buffered_polygons_filtered_sf <- sf::st_sfc(buffered_polygons_filtered, crs = sf::st_crs(seeded_points)) |> st_as_sf()
# Intersect the plots with the reef polygons, calculate overlap
intersections_filtered <- st_intersection(buffered_polygons_filtered_sf, st_union(geomorphic_seeding_filtered)) %>%
mutate(intersection_area = st_area(.)) %>%
mutate(percentage_overlap = (as.numeric(intersection_area)/10000) * 100) |>
mutate(overlap=ifelse(percentage_overlap==100, "white", "darkred"))
toc()
## 2.952 sec elapsed
tmap_mode("plot")
polyplots <- tm_basemap("Esri.WorldImagery") +
tm_shape(geomorphic_seeding_filtered, bbox=inset) +
tm_polygons(fill = "class",
fill.legend = tm_legend(title="Habitat", position = tm_pos_out("right", "center")),
fill.scale=tm_scale_categorical(values=habitat_pal_filtered),
lwd = 0.2, col="black", fill_alpha = 1) +
tm_shape(intersections_filtered,
bbox=inset) +
tm_polygons(fill=NA,
col="overlap",
lwd=0.5,
col_alpha=0.5,
col.scale=tm_scale(values=overlap)) +
tm_title("2) Seeded polygons", size=1)
polyplots
Use future.lapply to speed things up with larger
datasets (note: multisession below, use
multicoreon base R to obtain optimal results)
options(future.rng.onMisuse = "ignore")
library(future.apply)
tic()
# Set up parallel processing plan
plan(multisession, workers = availableCores())
# Parallel lapply sampling process using future_lapply
seeded_points_list <- future_lapply(future.seed=NULL, seq_len(nrow(geomorphic_seeding)), function(i) {
polygon <- geomorphic_seeding[i, ]
points <- st_sample(polygon, size = polygon$n, type = "random")
st_sf(class = polygon$class, geometry = points)
})
# Combine the list into a single sf object
seeded_points <- do.call(rbind, seeded_points_list) |> st_sf() |> st_transform(20353)
# Define the width and length of the rectangle
width <- 100 # Example width (in the same units as your CRS)
length <- 100 # Example length (in the same units as your CRS)
# Use future_lapply to create the buffered polygons
buffered_polygons_list <- future_lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
return(polygon)
})
# Combine the results into a single sf object
buffered_polygons_sf <- sf::st_sfc(buffered_polygons_list, crs = sf::st_crs(seeded_points)) |> st_as_sf()
# Intersect the plots with the reef polygons, calculate overlap
intersections <- st_intersection(buffered_polygons_sf, st_union(geomorphic_seeding)) %>%
mutate(intersection_area = st_area(.)) %>%
mutate(percentage_overlap = (as.numeric(intersection_area) / 10000) * 100)
# Reset the plan to sequential processing after completion (optional)
plan(sequential)
toc()
Determine optimum plot 1ha locations using polygons with randomly seeded centroids
rotate polygons by a random 5° compass bearing between 0-355°.
retain polygons that are 100% overlap with habitat type.
use future.lapply to multisession
library(future.apply)
library(sf)
library(spatialEco)
tic()
# Set up parallel processing plan
plan(multisession, workers = availableCores())
# Use future_lapply to create the buffered polygons with random rotations
buffered_polygons_list_rotated <- future_lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
# Convert to sf object for rotation
polygon_sf <- st_sf(geometry = st_sfc(polygon, crs = st_crs(seeded_points)))
# Select a random rotation angle from the sequence
angle <- sample(seq(0, 355, 5), 1)
# Rotate the polygon by the random angle using spatialEco::rotate.polygon
rotated_polygon <- spatialEco::rotate.polygon(polygon_sf, angle = angle, anchor = "center") |> mutate(id=i)
# Extract the geometry from the rotated sf object
return(rotated_polygon)
})
buffered_polygons_rotated <- do.call(rbind,buffered_polygons_list_rotated) |> st_set_crs(20353)
buffered_polygons_rotated_intersect <- st_intersection(buffered_polygons_rotated, st_union(geomorphic_seeding_filtered)) %>%
mutate(intersection_area = st_area(.)) |> as.data.frame() |> select(-geometry)
intersections_rotated <- buffered_polygons_rotated |>
left_join(buffered_polygons_rotated_intersect, by="id") |>
mutate(percentage_overlap = (as.numeric(intersection_area) / 10000) * 100)
# Reset the plan to sequential processing after completion (optional)
plan(sequential)
toc()
## 18.816 sec elapsed
inset <- st_bbox(c(xmin = 145.44, xmax = 145.455, ymin = -14.697, ymax = -14.692)) |>
st_bbox()
tmap_mode("plot")
rotatedplots <- tm_basemap("Esri.WorldImagery") +
tm_shape(geomorphic_seeding_filtered, bbox=inset) +
tm_polygons(fill = "class",
fill.legend = tm_legend(title="Habitat", position = tm_pos_out("right", "center")),
fill.scale=tm_scale_categorical(values=habitat_pal_filtered),
lwd = 0.2, col="black", fill_alpha = 1) +
tm_shape(intersections_rotated |> filter(!percentage_overlap>98),
bbox=inset) +
tm_polygons(fill=NA,
lwd=0.5,
col_alpha=0.5,
col="darkred") +
tm_shape(intersections_rotated |> filter(percentage_overlap>98),
bbox=inset) +
tm_polygons(fill=NA,
lwd=1,
col="white") +
tm_title("3) Rotated random polygons", size=1)
rotatedplots
Broadly similar n between rotated and parallel polygons:
combined <- rbind(intersections_filtered |> mutate(type="Parallel plots") |> as.data.frame() |> select(intersection_area, type),
buffered_polygons_rotated_intersect |> mutate(type="Rotated plots") |> as.data.frame() |> select(intersection_area, type)
)
combined <- combined |> mutate(percentage_overlap=as.numeric(intersection_area)/10000*100)
combined_data <- as.data.frame(combined)
hist_data <- ggplot_build(ggplot(combined_data) + geom_histogram(aes(x = percentage_overlap, fill = type)))$data[[1]]
scale_factor <- max(hist_data$count) / max(hist_data$density)
ggplot() +
theme_bw() +
facet_wrap(~type) +
stat_density(data = combined_data, aes(x = percentage_overlap, y = ..density.. * scale_factor, col = type),
alpha = 0.2, linewidth = 0.5, show.legend=FALSE, fill=NA) +
geom_histogram(data = combined_data, aes(x = percentage_overlap, y = ..count.., fill = type),
alpha = 0.2, color = "black", linewidth = 0.2, show.legend=FALSE) +
scale_y_continuous(
name = "Count",
sec.axis = sec_axis(~ . / scale_factor, name = "Density")
) +
scale_fill_manual(values=c("turquoise4", "coral4")) +
scale_color_manual(values=c("turquoise4", "coral4"))
No difference in n plots between methods, but: rotating plots improves fit in narrow irregular polygons: Comparing the outputs from the three approaches:
tmap_arrange(gridplots, polyplots, rotatedplots, ncol=1)
sf is very slow for large computations. Opt for
shapely or geopandas:
from shapely.geometry import Polygon
from shapely.affinity import rotate
import geopandas as gpd
from shapely.affinity import rotate